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1.  INTRODUCTION 
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In  a recent  study  [1]*  of  the  convergence  properties  of  finite  element 
methods  in  nonlinear  fluid  mechanics,  an  indirect  approach  was  taken.  A two- 
dimensional  example  with  a known  exact  solution  was  chosen  as  the  vehicle  for 
the  study,  and  various  mesh  refinements  were  tested  in  an  attempt  to  extract 
information  on  the  effect  of  the  local  Reynolds  number.  However,  more  direct 
approaches  are  usually  preferred.  In  this  study  one  such  direct  approach  is 
followed,  based  upon  the  spectral  decomposition  of  the  solution  operator. 

Spectral  decomposition  is  widely  employed  as  a solution  technique  for 
linear  structural  dynamics  problems  and  can  be  applied  readily  to  linear, 
transient  heat  transfer  analysis;  in  this  case,  the  extension  to  nonlinear 
problems  is  of  interest.  It  was  shown  in  [2]  that  spectral  techniques  were 
applicable  to  stiff**  systems  of  rate  equations,  while  recent  studies  [3,4] 
of  geometrically  and  materially  nonlinear  structural  dynamics  have  demon- 
strated the  increased  information  content  of  the  numerical  results.  The  use 
of  spectral  decomposition  in  nonlinear  problems  of  heat  and  mass  transfer 
would  be  expected  to  yield  equally  increased  flow  of  information  to  the 
analyst,  and  this  information  could  include  a quantitative  comparison  of  various 
solution  strategies,  meshes,  and  element  hierarchies. 

In  order  to  examine  the  use  of  spectral  techniques  in  nonlinear  heat 
and  mass  transfer,  a relatively  simple  framework  was  chosen,  based  upon  the 
description  of  advection  (convection)  and  diffusion  (viscosity).  Advection- 
diffusion  equations  describe  numerous  transport  phenomena  of  interest  to 
engineers  and  scientists.  These  descriptions  typically  have  the  form 

+ “it,'  ■ • ti*»  • <’•’> 


* Brackets  denote  references  listed  at  the  end  of  the  paper. 

**  Stiff  systems  are  those  with  a large  spread  in  the  eigenvalue  spectrum  of 
the  governing  matrix  operater. 


where  $ is  the  system  quantity  being  transported,  p is  the  density,  a the 
diffusion  parameter,  t the  time,  and  u.  the  Cartesian  components  of  advection 
velocity  in  the  spatial  coordinate  directions. 

For  realistic  situations  (1.1),  with  an  appropriate  set  of  initial  and 
boundary  conditions,  describes  a problem  that  is  too  complex  for  analytic 
treatment.  Therefore,  (1.1)  is  usually  solved  approximately  by  some  numerical 
method,  such  as  the  methods  of  finite  difference  or  finite  element.  When  a 
continuous  partial  differential  equation,  such  as  (1.1),  is  reduced  to  an 
"equivalent"  set  of  discrete  equations,  the  discrete  set  is  found  to  display 
characteristics  unlike  those  found  in  the  original  equation.  Such  charac- 
teristics include  both  temporal  and  spatial  oscillations,  artificial  dissi- 
pation, and  dispersion.  Since  this  behavior  arises  strictly  from  the  numerical 
approximation,  it  must  be  understood  and  controlled  in  order  for  the  approxi- 
mate methods  to  achieve  any  validity  and  utility. 

Finite  difference  approximations  of  the  advection-diffusion  equation 
have  been  investigated  widely  (see  [5],  e.g.)  and  continue  to  be  an  area  of 
research  interest.  Finite  element  approximations  to  (1.1)  are  of  more  recent 
vintage  [6],  but  are  receiving  a good  deal  of  attention  in  the  current  litera- 
ture (e.g.,  see  [7,8]).  However,  little  in  the  way  of  analysis  of  the  numeri- 
cal characteristics  is  available  [9].  In  the  following  sections  the  discretized 
equations  are  examined  through  the  techniques  of  spectral  decomposition,  as 
applied  to  the  one-dimensional  prototype  of  (1.1).  Although  the  example  is 
one-dimensional,  all  of  the  appropriate  characteristics  of  more  complex 
examples  are  present  for  study. 
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2.  ONE-DIMENSIONAL  TRANSPORT  EQUATIONS 


The  analysis  of  the  multi-dimensional  equation  given  in  (1.1)  is  too 
complex,  even  in  discrete  form;  thus,  several  one-dimensional  cases  will  be 
considered.  One  of  the  popular  nonlinear  forms  of  (1.1)  which  is  often 
analyzed  is  the  Burger's  equation 
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0 < a < oo 


(2.1) 


Equation  (2.1)  is  a model  for  certain  aspects  of  turbulence  and  shock  wave 
studies;  it  is  also,  in  some  sense,  a one-dimensional  analogue  to  the  Navier- 
Stokes  equation  (though  there  is  no  pressure  gradient  term  or  incompressi- 
bility constraint). 

A second  equation  that  has  received  some  study  is  the  linear  advection- 
diffusion  equation  (color  equation) 


3<j) 
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0 < a<  oo. 


(2.2) 


Equation  (2.2)  is  a model  for  the  one-dimensional  transport  of  some  intensive 
property  of  a system;  e.g.,  energy  or  species  concentration. 

Both  of  the  above  equations  will  be  treated  in  the  following  sections 
to  illuminate  both  the  similarities  and  differences  (e.g.,  linear  vs.  non- 
linear). Also,  both  steady  and  transient  cases  will  be  considered.  In 
studying  (2.1)  and  (2.2),  it  must  be  emphasized  that  conclusions  arrived 
at  do  not  necessarily  carry  over  to  the  multi -dimensional  case.  However, 
it  is  felt  that  the  study  of  (2.1)  and  (2.2)  will  give  some  insight  into 
appropriate  solution  methods  and  modeling  criteria  for  the  higher 
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dimensionality  problems. 

To  complete  the  specification  of  the  model  problems,  suitable  initial 
and  boundary  conditions  for  (2.1)  and  (2.2)  are  required.  For  both  cases 
the  spatial  domain  must  be  finite  (for  computational  convenience  and  accuracy) 
and  is  chosen  to  be  a unit  length  on  the  positive  x axis;  i.e. , 0 < x < 1. 

For  the  Burger's  equation,  the  following  boundary  and  initial  conditions  will 
be  treated: 


(Steady  State)  u17=af^  » 0<x<l,  (2.3) 

with  u(l ) = 0,  u(0)  = 2; 

(Transient)  §£  + “ H = « 0 . 0 < x < 1 . (2.4) 

with  u (1  ,t ) = 0,  u(0,t)  = 2, 
and  u(x,0)  = 0; 


For  the  advection-diffusion  equation,  the  following  boundary  and  initial 
conditions  will  be  treated: 


(Steady  State)  u ff  = « 0 , 0 < x < 1 , (2.5) 

with  4>(1 ) = 0 , 4>(0)  =2  ; 


and  u = u(x)  specified. 
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(Transient)  ft  + u If  = “ ’ 0<x<l,  (2.6) 

with  4>(1  ,t)  = 0,  4>(0,t)  = 1 , 

*(x,0)  =0  , 

and  u = u(x)  specified. 
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3.  THE  EIGENSPECTRUM  PROBLEM 

The  discretization  of  the  continuous  Burger's  equation  leads  to  the 
governing  matrix  equation  (see  Appendix  A) 

M • u + jjc  + C(u)j  • u = F , (3.1) 

where  M is  the  linear  translational  inertia  matrix;  K is  the  diffusion  matrix; 

ss  ss 

C is  the  advective  matrix;  F is  the  vector  of  nodal  point  forces,  due  to 

*5*  ■** 

initial  and  boundary  values;  and  u is  the  vector  of  nodal  point  velocities. 

The  matrix  M is  symmetric,  positive-definite;  K is  symmetric  and  at  least 
semi-definite;  and  C is  generally  unsymmetric.  Although  (3.1)  is  written  for 
Burger's  equation,  it  should  be  pointed  out  that  the  same  matrix  equation 
characterizes  transient,  nonlinear  Navier-Stokes  flow.  Also,  the  development 
that  follows  is  equally  valid  for  the  more  general  case. 

The  generalized  eigenvalue  problem  [10]  is 


(3.2) 


* * 
where  $.  are  the  right-hand  eigenvectors  and  the  eigenvalues.  The  matrix 

G is  the  sum  of  the  diffusion  (viscosity)  and  the  instantaneous  advection 

(convection)  matrices.  Superposed  asterisks  indicate  a complex  scalar  or 

vector. 

Similarly,  for  the  transposed  eigenvalue  problem. 


* 


as 


(3.3) 
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where  iJk  are  the  left-hand  eigenvectors. 

The  solution  vector,  u,  can  be  expanded,  in  terms  of  either  set  of 
eigenvectors,  as 


u(t)  = y>;(t),;(t) 

i=l 


(3.4) 


or 


u(T)  -^(t^t). 


(3.5) 


i=l 


where  the  complex  modal  coefficients  are  ai  and  , respectively.  The  explicit 
time  dependence  of  the  eigenvalues  and  eigenvectors  indicates  the  changing 
eigenspectrum  of  the  nonlinear  system. 

Othogonality  implies  that 


ifj*  • M • 4> 

~j  * ~i 


0 , i i j 


mi  , i=j 


(3.6) 


and 


(V1'1 


(3.7) 
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If  the  conjugate  of  a complex  scalar  or  vector  is  denoted  by  a superposed 
bar,  then  the  modal  coefficients  at  any  instant  can  be  found  from  the  current 
velocity  field  as 


* r Hi- 

« T ( t ) = • M • u(t)  -J- 


* _* 
m3mj 


(3.8) 


A similar  expression  could  be  obtained  for  B*(t)  if  the  expansion  (3.5)  is 
chosen. 

Using  the  orthogonality  relations  (3.6)  and  (3.7),  the  uncoupled 
(modal)  Burger's  equations  are 


* • * ★ ★ * 

"j  “j  + 9j  “j  = fj. 


(3.9) 


where 


fj  = *:T-  e • 


(3.10) 


For  the  case  where  the  forcing  function  can  be  adequately  described  by  a 
linear  representation  over  the  time  interval  At  = tn+1-tp,  the  exact  solution 
for  each  modal  coefficient  is  given  by  the  sum  of  three  terms: 


aj(tn+l)  = aj0  + ajl  * At  + ai?  * 


-X*At 


(3.11) 


where 
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(3.12) 


(3.13) 


and 


•k  k 

“J2  ' “j(tn> 


(3.14) 


k k k 

The  quantities  aj(tn)>  fjQ,  and  f ^ represent  the  initial  condition  for  the 
mode,  the  initial  value  of  the  modal  forcing  function,  and  the  incremental 
addition  to  the  modal  forcing  function,  respectively.  The  initial  condition 
is  found  from  the  solution  at  time  tn  through  projection  into  modal  space: 


* *T 

a At  ) = ipi  • M • u(t  ) • -j br  J 
J ~ ~ m.m. 

J J 


(3.15) 


the  forcing  function  can  be  similiarly  found: 


fj(x)  = * F(t),  tn  < t < tn+1 


(3.16) 
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* * 1 
f-1  = Ip- 
jl  ~J 


AF 


(3.19) 


For  the  case  of  a forcing  function  that  is  constant. 
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(3.20) 


CjCJ 


It  should  be  pointed  out  that  the  complex  frequencies,  X.;  eigenvectors, 
★ ★ J 

and  <j).;  and  other  modal  quantities  are  theoretically  changing  continuously 

vJ  "“J 

with  time,  since  the  problem  is  nonlinear.  This  development  is  based  upon 
small  departures  from  the  nonlinear  state  during  a particular  time  step,  how- 
ever; methods  for  extrapolating  the  eigenspectrum  could  be  studied,  but  such  a 
study  is  beyond  the  scope  of  this  investigation. 
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4.  ANALYTIC  SOLUTIONS  FOR  MODEL  EQUATIONS 


For  purposes  of  comparison  and  determination  of  solution  accuracy  it 
is  convenient  to  have  a closed  form  analytic  solution  to  the  test  problems. 
For  the  Burger's  equation  a variety  of  analytic  solutions  are  available  for 
both  steady  and  transient  problems,  and  for  a wide  variety  of  boundary  and 
initial  conditions.  These  solutions  have  been  cataloged  by  Benton  and 
Platzman  [11].  For  the  cases  considered  here,  the  analytic  solutions  may  be 
expressed  by: 

(Steady  State)  u(x)  = -A  tanh  (4.1) 


where  B = 1 

1 = tanh(^)  ; 


(Transient) 


u(x,t)  = 2 


(4.2) 


where 


F(x,t)  = 1 et-x 


erfc 


~x-2t 

2/r 


In  the  last  expression  for  F(x,t),  a change  of  scale  and  a shift  of  origin 
are  required  to  allow  u(x,t)  to  describe  the  specified  problem. 
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5.  NUMERICAL  EXAMPLES 


Analyses  of  Burger's  equation  were  carried  out  for  a wide  range  of  the 
diffusion  parameter,  a,  while  maintaining  identical  initial  and  boundary 
conditions  on  the  velocity  field.  Three  cases  were  selected  as  representative: 
a = 50,  where  diffusion  dominates;  a = 1,  where  diffusion  and  advection  are 
balanced;  and  a = .01,  where  advection  dominates.  For  each  value  of  a,  three 
solutions  were  obtained.  The  first,  obtained  from  the  closed-form  expression 
(4.2),  will  be  referred  to  as  the  "exact"  solution.  The  second  was  obtained 
from  an  implicit,  Crank-Nicolson  integration  of  the  discretized  equation  (3.1), 
and  will  be  termed  the  "direct  integration"  solution  (the  procedure  is  des- 
cribed in  Appendix  B).  The  third  was  the  result  of  carrying  out  the  eigen- 
spectrum  analysis  of  the  discretized  matrix  equation,  followed  by  "precise" 
eigenmode  integration,  at  each  time  step.  This  solution  will  be  called  the 
"eigenspectrum"  solution.  The  eigensystem  was  computed  using  an  algorithm 
developed  by  Moler  and  Stewart  [12]. 

The  major  findings  of  this  study  are  exhibited  in  Figures  1,2,  and  3. 

For  each  method  (exact  solution,  direct  integration,  mode  superposition),  four 
distinct  solutions  were  obtained,  with  order  of  magnitude  changes  in  the  time 
step.  Ten  steps  were  taken  for  each  of  the  four  solutions.  To  the  extent  that 
the  distinct  solutions  are  continuous  from  one  decade  to  the  next,  the  time 
step  can  be  interpreted  as  being  too  small  - namely,  unless  the  early  time 
information  has  some  intrinsic  value,  it  should  not  be  computed,  since  it  is 
not  needed  to  advance  the  solution  accurately.  We  have  chosen  to  plot  the  mid- 
side node  of  the  first  quadratic  element  - that  node  that  is  most  near  the 
initial  velocity  discontinuity  that  propagates  downstream.  This  midside  node 
sees  the  greatest  variation  in  velocity  at  the  earliest  time  and  should  pro- 
vide the  most  sensitive  measure  of  accuracy. 

For  a = 50  (Figure  1),  diffusion  is  the  dominant  transport  mechanism. 

-4  -3 

The  problem  is  solved  for  a time  step  of  10  from  0.  to  10  ; again,  for  a 

-3  -2  -2 

time  step  of  10  from  0.  to  10  ; again,  for  a time  step  of  10  from  0.  to 

10~^;  and,  finally,  for  a time  step  of  10”^  from  0.  to  1 . The  solution  for 

the  last  decade  is  plotted  on  a shifted  scale.  This  can  be  deduced  by 


Figure  2.  Comparisons  for  ct=1(L=3^. 


observing  that  the  exact  solution  is  always  continuous  from  decade  to  decade. 
Note  that,  for  the  first  two  decades,  both  the  direct  integration  and  modal 
superposition  methods  are  virtually  identical,  and  that  both  are  subject  to 
substantial  dispersion.  This  dispersion  (or  artificial  momentum  transport) 
is  due  to  the  lack  of  sufficient  spatial  definition  at  early  time.  The  direct 
integration  method  displays  a slight  discontinuity  between  second  and  third 
decades  (note  that  the  mode  superposition  remains  continuous),  with  a large 
discontinuity  at  the  beginning  of  the  final  decade.  The  mode  superposition 
method  has  a modest  discontinuity  between  the  third  and  fourth  decades. 

Based  upon  these  data,  the  mode  superposition  method  can  be  used  with  a time 
step  about  one  order  of  magnitude  larger  than  that  for  the  Crank-Nicolson 
method,  while  maintaining  the  same  accuracy,  when  diffusion  dominates. 

For  a = 1 (Figure  2),  a distance  equal  to  three  times  the  shock  thick- 
ness (the  shock  thickness  is  equal  to  a)  was  modeled  with  ten  quadratic 
elements.  The  same  sequence  of  time  steps  and  decades  as  that  used  for  the 
a = 50  studies  was  repeated  for  this  case  of  balanced  diffusion  and  advection. 
Again,  artificial  dispersion  dominates  the  first  two  decades.  The  mode 
superposition  method  is  slightly  less  accurate  for  the  first  two  decades; 
however,  between  the  second  and  third  decades,  the  direct  integration  solution 
suffers  a discontinuity,  while  the  mode  superposition  becomes  extremely 
accurate.  The  fourth  decade  is  plotted  on  a shifted  scale,  again,  and  both 
methods  become  less  accurate;  mode  superposition  is  the  most  rapidly  converg- 
ing of  the  two  methods  during  this  period.  As  before,  mode  superposition 
can  be  used  with  a time  step  about  one  order  of  magnitude  greater  than  that 
for  direct  integration,  even  here  where  diffusion  and  advection  are  both  of 
similar  importance. 

Finally,  for  a = .01  (Figure  3),  a distance  equal  to  five  shock  thick- 
nesses is  modeled  with  ten  quadratic  elements.  The  problem  is  solved  for 
-6  -5  -4  -3 

time  steps  of  10  ,10  ,10  , and  10  , with  each  decade  taking  ten  time 

steps.  For  such  a small  value  of  a the  advection  terms  dominate  the  trans- 
port of  momentum.  As  before,  and  seemingly  for  all  values  of  a,  artificial 
dispersion  causes  both  numerical  solutions  to  fall  well  below  the  exact 
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solution  during  the  first  two  decades,  with  the  direction  integration  solution 
slightly  less  dispersive.  Both  numerical  solutions  approach  the  exact  solution 
in  the  third  decade  - the  mode  superposition  results  become  asymptotic  and 
essentially  exact,  while  the  Crank-Nicol son  results  drift  above  the  exact 
solution.  This  latter  trend  is  considerably  exaggerated  in  the  fourth  decade, 
where  the  superior  convergence  properties  of  mode  superposition  are  easily 
observed.  Again,  about  an  order  of  magnitude  increase  in  the  time  step  can 
be  tolerated  with  the  mode  superposition  method  over  that  required  for  the 
Crank-Nicol son  approach. 

Spatial  convergency  can  be  studied  in  several  ways.  Two  methods  are 
described  here.  First,  the  analyst  can  experiment  with  a variety  of  dis- 
cretizations, while  noting  the  point  convergence  of  the  velocity  field.  For 
a = .01,  the  steady-state  numerical  solution  (mode  superposition  and  direct 
integration  gave  identical  results)  is  compared  to  the  exact  solution,  as  a 
function  of  the  number  of  elements  in  the  mesh  (Table  1).  The  elements  are  all 
quadratic.  As  can  be  seen,  even  the  points  within  the  shock  thickness  are 
adequately  treated  by  the  five-element  model.  Another,  and  more  informative, 
approach  is  to  examine  the  eigenvalue  spectrum  as  a function  of  mesh  size. 

For  a = 50,  eigenvalue  convergence  is  demonstrated  for  a variety  of  element 
sizes  (Table  2).  These  eigenvalues  were  calculated  at  an  early  time  in  the 
solution  but,  unlike  many  structural  problems,  the  spectrum  is  altered  only 
slightly  by  system  non-linearity.  For  this  reason,  iterative  methods  for 
computing  and  recomputing  the  eigenvalue  spectrum  should  be  extremely  economical. 
This  same  statement  may  not  be  valid  for  two-dimensional  flows  with  more  complex 
structure  (e.g.,  recirculation  regions,  separation  and  reattachment,  etc.). 

The  eigenvalues  and  eigenvectors  are  easily  estimated  (see  Appendix  C, 
Equation  (C.6)).  Table  3 shows  the  comparison  between  the  estimated  values  of 
the  characteristic  times  and  the  computed  values  for  a ten-element  model  for 
each  of  the  three  a regimes.  Only  the  first  twelve  modes  (of  the  nineteen 
non-rigid-body  modes  in  the  model)  are  estimated.  In  all  cases,  the  esti- 
mated values  are  almost  identical  with  the  computed  values  for  the  first  ten 
modes;  however,  the  caluclated  mode  shapes  begin  to  deteriorate  rapidly  after 
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TABLE  2.  Eigenvalue  convergence  (a=50). 
Number  of  elements,  N. 


Eigenvalue,  X 

N=2 

N=5 

N=10 

N=20 

1 

.2218 

.2176 

.2174 

.2174 

2 

.8200 

.8134 

.8099 

.8096 

3 

2.5648 

1.8269 

1.7988 

1.7967 

4 

3.3167 

3.1896 

3.1790 

5 

5.0201 

4.9943 

4.9575 

6 

8.6743 

7.2358 

7.1340 

7 

13.1623 

9.9503 

9.7114 

Highest 

2.5648 

26.4666 

116.1470 

476.0578 
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the  tenth  mode,  often  displaying  an  inconsistent  number  of  changes  in  sign 
(crossings)  in  the  modal  shape.  As  an  example,  the  first  mode  right-  and  left- 
hand  eigenvectors,  for  a = .01,  normalized  to  a maximum  amplitude  of  unity, 
are  depicted  in  Figure  4.  Note  the  rapid  rise  of  the  vectors  within  the 
shock  thickness  at  each  end  of  the  mesh.  The  fact  that  the  two  vectors  are 
nearly  mirror  images  of  each  other,  for  this  uniform  mesh  of  ten  quadratic 
elements,  may  be  of  importance  in  explaining  the  effects  of  varying  mesh 
size,  as  will  be  discussed  later. 

In  order  to  examine  the  cause  for  deterioration  of  the  accuracy  of  the 
eigenspectrum  after  about  ten  modes,  the  tenth  mode  right-  and  left-hand 
eigenvectors  are  plotted  in  Figure  5.  Note  that,  for  this  mode,  the  corner 
and  midside  nodes  for  the  quadratic  element  almost  perfectly  coincide  with 
relative  maxima,  relative  minima,  and  null  points  for  the  modal  shape. 

Higher  modes  than  this  exceed  the  capacity  of  the  quadratic  element  to  fit  a 
second-order  polynomial  to  known  information  concerning  these  critical  points. 
Higher-order  elements  would  be  expected  to  continue  to  track  the  eigenspec- 
trum until  their  modeling  capacity  was  exceeded.  This  provides  us  with  a 
method  for  estimating  the  number  of  retained  modes  in  an  analysis,  however. 

In  one  dimension,  if  the  element  length  is  2Ax  (ax  being  the  distance 
between  a corner  node  and  a midside  node),  the  n = L/2ax. 

It  is  not  always  feasible  to  use  a uniform  mesh  for  production  analysis. 
Therefore,  a modest  examination  of  one  other  modeling  feature  was  undertaken. 
The  graded  mesh  seemed  to  exhibit  the  undesirable  property  of  causing  the 
eigenspectrum  to  become  indistinct,  bunching  the  eigenvalues  into  pairs. 

Under  certain  conditions,  complex  eigenvalue/eigenvector  pairs  were  pro- 
duced. Occasionally,  these  complex  pairs  appeared  early  in  the  solution, 
disappearing  as  the  solution  approached  the  steady  state.  In  all  cases  to 
date,  when  the  mesh  was  carefully  prepared,  these  characteristics  were  not 
observed.  Although  further  study  is  needed,  one  can  speculate  that  several 
events  are  coalescing  in  order  to  produce  this  behavior:  (1)  a is  small, 
indicating  that  advection  is  the  important  transport  mechanism;  (2)  the 
velocity  is  relatively  uniform,  leading  to  a skew-symmetric  form  for  the 
advection  matrix;  and  (3)  grading  of  the  mesh  can  lead  to  difficulty  in 
calculating  the  "mirror  image"  left-  and  right-hand  eigenvectors. 
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DISTANCE,  X 

Figure  4.  First  mode  eigenvectors  for  a=0.01. 


6.  CONCLUSIONS 


\ ' 

A spectral  decomposition  method  based  upon  finite  element  modeling  has 
been  compared  to  a Crank-Nicol son  direct  integration  solution  scheme  and  the 
exact  solution  for  the  one-dimensional,  nonlinear  system  defined  by  Burger’s 
equation.  Results  from  this  study  are  applicable  to  both  fluid  mechanics 
and  combined  conduction-convection  heat  transfer.  The  parameter  a,  v/hich 
governs  the  importance  of  diffusive  transport,  was  varied  over  a sufficiently 
wide  range  such  that  comments  on  the  comparisons  are  general. 

The  mode  superposition  method  proved  to  be  very  attractive,  in  com- 
parison to  the  second-order  accurate  Crank-Nicol son  approach,  generally 
allowing  an  order  of  magnitude  larger  time  step  for  equivalent  convergence 
to  the  exact  solution.  The  modal  shapes  themselves  tend  to  provide  useful 
information  about  the  ability  of  a given  mesh  to  produce  accurate  results, 
much  in  the  same  way  that  modal  information  is  used  in  nonlinear  structural 
dynamics.  For  this  class  of  problems,  in  contrast  to  structural  dynamics, 
system  nonlinearities  did  not  manifest  themselves  in  dramatic  changes  in 
the  eigenspectrum. 
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APPENDIX  A 


FINITE  ELEMENT  MODELS  FOR  BURGER’S  EQUATION 
Consider  now  the  transient  form  of  the  Burger's  equation 


9u 

at 


+ 


ax 


„a2u 

*djr 


(A.  1 ) 


It  is  desired  to  perform  a spatial  discretization  on  the  above  equation 
using  the  finite  element  method  in  conjunction  with  a Galerk.in  procedure. 

Let  the  dependent  variable  be  approximated  by 


u(x,t)  = 2$i(x)  • u.(t)  (A. 2) 

i 


or 


u(x,t)  = *T(x)  • u(t) 


(A. 3) 


Substituting  (A. 3)  into  (A.l)  yields 


at  au  A j ar>T  „ a2$T  D 

0 ~ + 0 u ~ u - a ~ u = R , 

at  ax  ax 


(A. 4) 
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where  R is  the  error  or  residual  resulting  from  the  approximation  in  (A. 3). 
The  Galerkin  procedure  makes  the  residual  zero  in  a weighted  sense  by 
forcing  R to  be  orthogonal  to  the  space  spanned  by  $>.  Thus, 


[ 

■ 


(A. 5) 


Equation  (A. 5)  may  be  rearranged  by  performing  an  integration  by  parts  on  the 
last  term,  i.e. 


/ 


3V 

- a<I>  ^5-  udV 


/3<3>T  f 3$  3?>T 

• at  sr  + J ir 


udV. 


Therefore,  (A. 5)  can  be  expressed  as 


>'«  ?•  \F--l 

- \/  J 1/ 


dV 


3?  3$ 

rt  9$T  \ 

+ 

/ <*-■=■  r=-  dV 

u = 

u ) dA 

J 3X  3X 

J ~\  3x  ~ / 

(A. 6) 
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In  Equation  (A. 6)  the  right  hand  side  represents  the  natural  boundary  condi- 
tion (i.e.  f-u-)  for  the  differential  equation. 

d X 

The  discrete  approximation  to  (A.l)  as  given  in  (A. 6)  can  be  completed 
once  the  basis  or  shape  functions  $(x)  are  specified.  For  the  present  case, 
the  basis  functions  to  be  considered  are  linear  and  quadratic  in  the  spatial 
coordinate.  Using  an  isoparametric  formulation  let 


Linear:  $(s)  j]/2(]+s)J 


l/2(s2-s) 

Quadratic:  4>(s)  =<  (1  -s2 ) 

(l/2(s2+s) 


(A. 7) 


(A. 8) 


where  s is  the  normalized  spatial  variable.  The  x coordinate  is  related 
to  s through 


x ( s ) = 1/2(1 -s)x  + l/2(l+s)x 
1 2 


(A. 9) 


where  x and  x are  the  coordinates  of  the  ends  of  the  element. 

1 2 

To  compute  the  previously  defined  integrals  of  the  shape  function 
derivatives  the  following  relations  are  required 


f|.  9X  ^ 

3s  as  3X 


or 


' ■ * « 1 So 

9x  J 3s 

ax  (x  -x  ) 

where  J = Jacobian  = ^ — from  (A. 9).  Thus, 

3$  2 3<J> 

3x  Ax  ¥s 

Also,  for  this  one-dimensional  case, 

I 

dV  = dx  ' 1 = ff  ds  * 1 = T-  ds 


and 


(A. 10) 


dV  = ds  • 1 . 


(A. 11) 


Linear  Elements 

Consider  first  the  case  of  linear  basis  functions,  i.e. 


*(s) 


\ l/2(l-s)  | 

1 1/2(1 +s ) j 
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Then  using  the  definitions  in  (A. 10)  and  (A. 11)  and  the  integrals  defined  in 
(A. 6)  allows  the  following  computations 


_r  r 


or,  symbolically. 


M u + C(u)  u + K u = F. 


(A. 14) 


Quadratic  Elements 

A development  similar  to  the  one  for  linear  elements  may  be  carried  out 
for  the  quadratic  basis  functions 


{l/2(s2-s)  ) 
(1-s2)  J 
l/2(s2+s) j 


The  element  matrices  are  then  given  by: 
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as  Ax  2 


(1/30) 


(-10u!-6u2+u3)  (12ui+  8u2) 
(-6ui-16u2+2u3)  (8ui-8u3) 


(-2ui-2u2-u3) 

(-2ui+162+6u3) 


(ui+2u2+2u3 ) (-8u2-12u3)  (-Ui+6u2+1 0u3 ) 


= C(u)  ; 
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T-rr 


and 
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K . (A. 1 5) 


This  results  in  a discrete  equation  for  each  element  of  the  following  form 


4/15  2/15  -1/15 
2/15  16/15  2/15 
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( -1  Ou  i -6u2"*’U  3 ) (12ui+8u2)(-2ui-2u2-u3) 

( -6u i -1 6u2+2u  3 ) (8u i -8u  3 ) ( -2u ! +1 6u  2+6u  3 ) 
(ui+2u2+2u3 ) (-8u2-12u3 ) (-Ui+6u2+1 0u3 ) 


(A. 16) 


or,  symbolically. 


M u + C(u)  • u + K u = F 


(A. 17) 
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appendix  b.  solution  strategy 

The  governing  system  that  defines  the  Crank-Nicol son  method 
of  the  matrix  Equation  (3.1),  evaluated  at  two  successive  times,  t 
Thus, 


K + 


• “n  + ? 


~n+l 


if  +£<W 


Hn+1 


= ~n  + !n+l 


If  a backward  difference  expression  is  used  to  eliminate  u and 

. ~n+l 

difference  expression  to  eliminate  u , 


~n+l 


_ Vl 


Un 


At 


= “n 


then  the  usual  Crank-Nicolson  method  is  derived, 


™ * Vl  + At 


P £<un+1) 


Vl  = At  En 


+ At  En+1  + 2,~1  ' % ' At 


K + C(un) 


~n 


is  the  sum 

, and  Vi- 


le.]) 


forward 


(B.2) 


(B.  3) 
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It  is  recognized  that  the  trapezoidal  rule  is  characteristically  equiva- 
lent to  the  Crank-Nicolson  method,  although  it  leads  to  a slightly  different 
form  for  Equation  (B.3).  Then, 


~n+l 


$(» 


n+1 


“n> 


(B.4) 


so  that  the  governing  equation  becomes 


K + 


?^n+l ) + At  5? 


-n+1 


-n+1 


+ M 

At  ~ 


*?n 


+ H 


(B.5) 


This  form  is  somewhat  more  convenient  in  terms  of  storage  for  the  forcing 
terms,  while  suffering  from  the  disadvantage  that  a hereditary  derivative 
appears  as  an  initial  condition. 

Both  methods  are  known  to  be  unconditionally  stable,  with  truncation 
error  of  order  (At)3.  Since  storage  was  not  considered  a problem  for  these 
one-dimensional  examples,  the  conventional  Crank-Nicolson  method  was  adopted. 

The  only  remaining  question  is  the  treatment  of  the  nonlinear  term 
C(un+1).  For  simplicity,  this  study  has  also  adopted  quasi-linearization: 


) 


(B.  6) 


It  is  recognized  that  this  approach  might  introduce  significant  errors  if 
the  nonlinearities  are  substantial.  In  such  a case,  a Newton-Raphson  or 
modified  Newton-Raphson  might  be  required  at  each  time  step. 
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APPENDIX  C.  EIGENVALUE/EIGENVECTOR  ESTIMATES 

In  order  to  arrive  at  estimates  of  the  eigenvalues  and  eigenvectors  of 
the  nonlinear  system  (2.1),  consider  the  linearized  version 


3u  . ?r  9u  _ a2u 

at  u 9x  ax7 


(C.l) 


with  u prescribed  on  the  boundaries  x=0  and  x=l , as  well  as  at  the  initial 
time  tQ.  Looking  for  solutions  of  the  form  e^<j>(x),  we  are  led  to  the 
eigenvalue  problem 


<j>"  - 2C<J>’  = A <p,  <p(0)  = <j>(l)  = 0. 


(C.2) 


In  the  pure  diffusion  case,  C=0,  and  the  eigenvalues  and  eigenvectors  are 
given  by 


Xn  = ~nZ^2  * ^ = sin(nirx) 


(C.  3) 


As  C increases,  corresponding  to  a decrease  in  a=l/2C,  the  skew-sym- 
metric advection  term  becomes  important  and  the  modes  change  form.  The 
eigenvalues  remain  real,  given  by 


A 


n 


-n2ir2-C2 


(C.4) 
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... 


while  the  eigenfunctions  are  given  by 


<j>n(x)  = e+Cxsin(n7rx) 


(C.5) 


The  plus  and  minus  signs  are  indicative  of  the  redefinition  of  orthogonality 
required  for  the  advection-diffusion  system,  with  the  left-hand  and  right- 
hand  eigenvectors  corresponding  to  functions  concentrated  at  opposite  bound- 
aries. 

For  our  problem,  it  is  convenient  to  write  Equation  (C.4)  in  the  form 


A_  = 


/ u2  n2n2\ 

-“(SF  + T7-)  • 


(C.6) 


where  a is  the  diffusion  parameter 
elements.  It  should  be  noted  that, 
stant  (for  our  problem,  uQ(x)=2  for 
ceeds,  however,  the  velocity  within 


and  L is  the  length  that  is  divided  into 
initially,  the  velocity  is  virtually  con- 
x>0  and  uQ(0)  =0).  As  the  analysis  pro- 
the  shock  thickness  varies. 
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